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ABSTRACT 


This thesis studies the fuel optimal periodic reboost profile required to maintain a 
spacecraft experiencing drag in low-earth-orbit (LEO). Recent advances in 
computational optimal control theory are employed, along with a Legendre-Gauss- 
Lobatto Pseudospectral collocation code developed at the Naval Postgraduate School, to 
solve the problem. Solutions obtained by this method are compared against a previous 
study. Key issues were checking the optimality of the solutions by way of the necessary 
conditions and the behavior of the solution to changes in the thruster size. The results 
confirmed Jensen's findings of propellant savings of one to five percent when compared 
against a middle altitude Forced Keplerian Trajectory (FKT). Larger savings are 
predicted if compared against a finite-bum Hohmann transfer with drag. The costates 
estimates compared favorably against necessary conditions of Pontryagin’s Minimum 
Principle. Analysis of the switching function yielded periods of thrust-modulated arcs. 
The optimal thrust profile appears to be a thrust-modulated bum to raise the orbit 
followed by an orbital decay and a terminating thrust-modulated arc. For a sufficiently 
low thrust-control authority, the switching structure includes a maximum thrust arc. 
Indirect optimization techniques to confirm these findings were unsuccessful. 
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I. 


INTRODUCTION 


Orbiting spacecraft experience a number of orbit perturbations, which usually 
require thruster firings to correct. These thruster firings use propellant, which is a non- 
replaceable resource in the spacecraft and often limits the mission duration. This thesis 
presents a numerical study of an optimal periodic thrusting method for low earth orbits 
for which drag is the primary orbital perturbation. 

Conventional thinking holds that the Hohmann transfer is the minimum energy 
transfer method, and hence, optimal. While this may be true for exoatmospheric orbits 
and ideal thrusters in which the impulse is applied instantaneously, it is not necessarily 
true for spacecraft in low earth orbit (LEO) with finite-bum thrusters. In fact, Ross and 
Alffiend have shown that there exists an orbit transfer method that is more efficient than a 
Hohmann transfer [Ref. 1]. Ross [Ref. 2] also showed that optimal endoatmospheric 
maneuvers generally contain “singular thrust arcs”. To quantitatively determine the 
optimal orbital maintenance maneuver, Jensen [Ref. 3] numerically investigated the 
problem based on algorithms developed by Fahroo and Ross [Ref. 4]. Part of this thesis 
is a follow-on to that analysis and seeks to confirm those findings using different 
numerical tools developed at the Naval Postgraduate School (NPS). 

An optimal orbital control methodology has potential to save thousands of dollars 
in launch costs and/or increase mission durations. This is particularly important 
considering the large number of spacecraft and constellations of spacecraft planned for 
low earth orbit. The propellant savings may be used directly to reduce spacecraft launch 
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mass at a savings of approximately $10,000 per kilogram [Ref. 3], or, the same amount of 
propellant may be launched, but the mission duration extended. Either of these are 
significant enough benefits to pursue this investigation. 

The optimal control analysis starts with the equations of motion and a cost 
function. The equations of motion describe the physical system. The cost function 
describes the amount of propellant consumed. The objective is to minimize the cost 
function consistent with the physical limitations of the system. 

Chapter II of this thesis contains a description of key concepts and methodologies 
along with the problem formulation for the optimal control problem. It first describes 
two different orbit transfer methods that are used for comparison to the optimal control 
method derived later. Optimal control theory is described along with a spectral 
collocation method used to discretize the problem for numerical analysis. Finally, the 
specific problem to be solved is formulated. 

Chapter III contains analysis performed on the use of linear versus non-linear 
equations in formulating orbital problems. This chapter studies and compares the use of 
Hill's linear equations of relative motion and the more typical nonlinear equations of 
motion. [Ref. 5] A well-known problem from Biyson and Ho is solved to confirm the 
solutions for the nonlinear equations. [Ref. 6] The linear equations are investigated to 
see if they provide a suitable replacement for the nonlinear equations. If so, this may 
benefit the numerical analysis by simplifying the equations of motion and possibly 
reducing computational time. 
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Chapter IV contains results from analysis using the direct method of optimization. 
The numerical analysis methodology is described along with the resulting optimal states, 
costates, and costs. The resultant states and controls are then compared against the 
necessary conditions described in Chapter n. Issues encountered during the numerical 
analysis conclude this chapter. 

Chapter V discusses the attempted implementation of an indirect method of 
optimization. A converged solutions was not obtained, so the “best” answer obtained is 
discussed along with the numerical analysis issues encountered. 

Chapter VI briefly looks at the thrust profiles. Some profiles encountered in this 
analysis were bang-bang while others followed singular arcs. This chapter relates those 
results to the switching function, which is obtained from the optimality conditions 
described in Chapter II. 

The thesis ends with a summary of major findings and conclusions. 
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II. OPTIMIZATION THEORY AND PROBLEM FORMULATION 


A. ORBITAL MAINTENANCE METHODS 

Spacecraft orbits are perturbed by a number of forces and the magnitude of these 
perturbations depend upon the specific orbit and spacecraft. This study looks at the 
impact of drag on low earth orbiting spacecraft and coplanar orbital transfer methods to 
counter the orbit decay caused by drag. 

Three different orbital maintenance methods are described here. These include 
the Hohmann Transfer, the forced Keplerian trajectory (FKT), and the periodic reboost. 
The Hohmann transfer boosts a spacecraft from one circular orbit to a different circular 
orbit using two boosts or thrustings. The FKT applies enough thrust to counter the drag 
so that thrust equals drag continuously. The periodic reboost does multiple bums to 
maintain boundary conditions and the number of bums are determined by a switching 
function. Both the Hohmann and periodic reboost method rely on boosting to a higher 
altitude and slowly decaying back to the original altitude at which time the maneuvers are 
repeated. 

The Hohmann transfer has long been considered the minimum energy or most 
efficient transfer method [Ref. 5]. It transfers a spacecraft between two orbits by using 
two tangential thrusts as shown in Figure II-1. The first bum, Av a , places the spacecraft 
into an elliptical orbit and the second bum, Avb, circularizes the orbit at the final altitude. 
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Figure II-1 Hohmann Transfer 


The FKT uses a drag cancellation process in which the thrust is continuously 
throttled to counter the force of drag. This requires the thruster to operate continuously at 
different thrust levels. Drag is a function of the local density (p), the orbital velocity (v), 
the area of the spacecraft in the direction of motion (A), and the coefficient of drag (Cd) 
(Drag=CDApv 2 /2). Since density varies during an orbit, the thrust level must be variable 
to exactly counter the varying drag force. 

The FKT varies with altitude and will be called either a low, mid, or high FKT. A 
low-FKT is an FKT performed at the initial altitude from which the Hohmann and 
periodic reboost transfers begin. The mid-FKT occurs at the altitude midway between the 
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initial and the highest orbit obtained by the Hohmann or periodic reboost. The high-FKT 
occurs at the highest altitude obtained by a comparable Hohmann or periodic reboost. 

The ideal Hohmann is closely approximated by the mid-FKT [Ref 1,3]- The low- 
FKT uses more fuel than the Hohmann while the high-FKT uses less fuel. Figure II-2 is 
from Jensen and gives a comparison of the propellant usage by each type [Ref. 3]. 

If there are no state constraints, it has been shown that the FKT is not the fuel 
optimal solution [Ref. 1]. Since the Hohmann transfer does not do better than the mid- 
FKT, the Hohmann reboost cannot be the fuel optimal solution either. This thesis 
attempts to identify an optimal periodic transfer that is more efficient. 


Propellant Comparison 



Figure II-2 Propellant Comparison from Jensen [from Ref. 3] 
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B. OPTIMIZATION THEORY 

The optimization problem is generally formulated by a system of state equations 
and a cost function. These equations are functions of the states (x), controls (u), and time 
(t). The problem is usually stated in the following manner. [Ref. 7] Given a dynamical 
system given by 

x = f(x,u,t) 


where boldface indicates vectors, determine the optimal control history, u*, which 
transfers the state of the system from its initial conditions to a final target, jp(Xf,tf) while 
minimizing the performance index (or cost), J, given by: 

J[u] = M(Xf,tf)+ j* f L(x,u,t)dt [Ref. 8] 

Pontiyagin's Minimum Principle provides the necessary conditions for optimality 
[Ref. 8]. The Hamiltonian, is constructed from the cost function (L) with the introduction 
of costates (X) given by: 

H = L(x,u,t) + X T f(u,x,t) 

The costates satisfy the following differential equation: 



The Minimum Principle states that the optimal control, u *, minimizes H at every 
point on the trajectory. For example, if there are no constraints on the controller then we 


must have 



and 


a 2 H 

du 2 


>0 


(n-3) 


The final conditions on the costates are obtained from the transversality equations 


[Ref. 8] 


M< f )A 

OX f 


fdv } 1 

K dX f ; 


(n-4) 




5 M (tf) , .,T dvttf) 


dtf 


■ +V 


at f 


= o 


(n-5) 


where \j/(x(tf),tf) = 0 defines the target states. 


C. PROBLEM FORMULATION 


The basic problem studied in Chapters IV through VI is a constrained 
optimization problem with both equality and inequality constraints. The cost function is 
of the Lagrangian form in which the cost is an integral in time. The five equations of 
motion detailed below in section 1 are the state equations of the form x = f where x is a 
vector containing the first order state equations. State constraints h(x) and control 
constraints g(u) also exist and may be equality (= 0) or inequality constraints (< 0). 

1. Normalized Equations of Motion 

The orbital equations of motion are first order ordinary differential equations. The 
equations are written for a coplanar low earth orbit in which drag (D) has a significant 
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effect on the orbit. A summary of the normalization process performed in Reference 3 is 
provided here. 

The first order equations of motion are based on the geometry in Figure II-3. 



The five equations of motion contain five states; radius (r), velocity (v), flight path 
angle (y), mass (m), and a reference angle (0), and two controls; thrust (T) and thrust 
angle (s) and are given by 

r = vsin(y) (H-6) 


. Tcos(s)-D . . . 

v=-—-gsin(y) 

m 


(D-7) 


r = 


/V ^ 

— g 


cos(y) T sin(^) 


(D-8) 


mv 


m = - 


-T 


(H-9) 


9 = 


v cos(y) 


(n-io) 
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When thrust (T) is normalized to a reference drag force, it allows a more intuitive 
interpretation of the results. A normalized thrust of one means the thrust equals the 
reference drag. A normalized thrust of 5 means the thrust is 5 times the force of the 
reference drag. 

The reference drag force is defined by the basic aerodynamic equation of drag 
given by 

D re f = ~Pref^-'D-^ V (P - ll) 


In this equation, the density p re f and the orbital velocity, v, are both at the 
reference altitude. The coefficient of drag, Cd, and area. A, are both physical 
characteristics of the spacecraft. 

The ballistic coefficient, B, of the spacecraft is also used to simplify the 
normalized equations. The ballistic coefficient is a function of the spacecraft's mass, 
coefficient of drag, and area (in the direction of velocity) and is given by 

B = -^- (n- 12 ) 

C d A 


If the ballistic coefficient is nondimensionalized by 

B- , B . 

j r refPref 

l 2 J 

then the equations of motion normalize to 
r = v • sin(/) 


v = -g • sin(y) + 


T • cos(^) - D 
mB 


cn-13) 


(n-14) 

(n-is) 


li 


(n-i6) 


y = 


7 y 2 J\ 
— “g 


cos(y) 


-T 
m =—= 

v e B 


— v 

0 = — • cos(y) 
r 


The state constraints for (r,v,y,m,0) are 

r (0) - r(tf) = 0 


-7t n 

-< y < — 

2 2 


T • sin(^) 
m-v • B 


v (0) - v (t f ) = 0 


(n-17) 

(11-18) 


Y (0) - Y (tf) = 0 


0 < m < 1 (normalized mass) 

0 <6 


The control constraints for (T,s) are 

0 < T < 5 (normalized thrust) 

-n <£<7C 

A generic spacecraft that experiences major disturbances due to drag is used in 
this study and is the same one used by Jensen [Ref. 3]. The specifications of this 
spacecraft, along with the normalized value where appropriate, are given in Table II-1. 
As a comparison, the normalized ballistic coefficients for some real spacecraft are 
included in Table 31-2. 
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Table II-1 Generic Spacecraft Characteristics from Jensen [from Ref. 3] 



Physical Units 

Normalized Units 

Area 

500 m 2 


Initial Mass 


1 

Maximum Thrust 

3.5 N 

5.0 

Initial Radius 

6678.15 km 

1 

Coefficient of Drag 

2.35 


Ballistic Coefficient 

2.55 kg/m 2 

40890 

Initial Orbit Radius 

6678.15 km 

1 


Table II-2. Ballistic Coefficients from Jensen [from Ref. 3] 


Generic Spacecraft 

4.089*10 4 

ISS-DACT 6 

1.26*10* 

Space Telescope 

4.72*10* 

Landsat-1 

4 04*10* 

Echo-1 

824 * 10 * 


2. Cost Function 

The objective is to minimize the amount of propellant required over a given 
period to maintain the orbit at or above the desired orbit radius. Therefore the cost is 
related to the change in mass divided by the change in time. The periodic cost function, 
J p , is given by 

m(0)-m(t f ) (n-l 9) 

tf 

where t f is the final time of the control period. Since thrust can be written as: 

T = -mv e (11-20) 

The cost function is rewritten as 

at (n-2i) 

V tr J V 
l f o v e 
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In terms of the non-dimensional (or normalized) variables, the periodic cost 


function becomes: 


J 


p 



(n-22) 


As discussed earlier in this thesis, the Hohmann and FKT orbit maintenance 
methods are not the most fuel efficient maneuvers. The cost function for an FKT 
trajectory is similar to the periodic cost function but with the additional refinement that 
the thrust equals the drag. So, the normalized thrust (normalized to drag) is equal to one. 
The FKT cost function then becomes: 


j fkt J=“Er dt CD-23) 

tf o v e B 

A ratio of these two cost functions provides an immediate indication of the 
performance of the periodic optimal control problem. If the ratio of the periodic cost to 
the FKT cost is less than one, then the periodic cost is more efficient than a low-FKT. 
The cost function then becomes. 


J = 



(D-24) 


D. NUMERICAL DISCRETIZATION 

A Legendre Pseudospectral method [Ref. 4] is used to formulate the periodic 

reboost problem for numerical analysis. This method uses polynomial approximations for 
the state and control functions and evaluates them at the Legendre-Gauss-Lobatto (LGL) 
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points. Discretizing the continuous problem into a finite dimensional nonlinear 
programming (NLP) formulation is necessary for numerical analysis. 

The LGL points lie within the interval [-1,1]. A transformation [Ref. 4] is used to 
change the cost, state, and final conditions from the interval [0, tf] to the interval [-1,1] 
resulting in the following cost and state equations 



(n-25) 


(H-26) 


or 




\ l f ) 


x(t) = f(x(t),u(t)) 


(H-27) 


This thesis addresses a periodic problem with periodic boundary conditions 
x(0)=x(tf). These boundary conditions become: 


x(-l) = x(l) (n-28) 

The state and control variables are approximated by N* order Lagrange 
interpolating polynomials on the interval [-1,1]. It can be shown [Ref. 4] that the cost 
function in Equation IT-25 can be rewritten in the following form. 

1 N 

J= oX w k- T k (D-29) 

z k=0 
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where w k are the LGL weights and T k are the values of thrust at each LGL point for 
k=0,...,N. 

The cost is now a discretized form of the integral in Equation 11-25. The state 
dynamics may be discretized as 

Ak = ~f( a k’^k) —c k 
N 

where c k = Z Z) kj a j for k - N, D kj are the elements of the differentiation 
j=0 

matrix, w k are the LGL weights, and a k , and b k are the values of the states and controls at 
t k , respectively. 

Similarly the system constraints can be approximated in the same manner. 

B k =g( a k,t>k)^0 ( n - 3 °) 

As can be seen from the above equation, this method of discretization retains 
much of the structure of the continuous problem and allows for easy numerical analysis. 
The code that implements this method is known as DIDO and was developed by 
Professors Fahroo and Ross at the Naval Postgraduate School (NPS). A front-end 
graphical user interface (GUI) was developed by Hall [Ref. 9] as part of his M.S. Thesis. 
In this thesis, both the GUI and non-GUI versions of DIDO were employed to simulate all 
the direct solutions. 
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III. LINEAR AND NONLINEAR ORBITAL EQUATIONS OF MOTION 


A. LINEAR ORBITAL EQUATIONS - HILL'S EQUATIONS 

1. Hill's Equations of Motion 

Hill's equations [Ref. 5] are useful in describing the relative motion between two 
close-orbiting satellites. The geometry is provided in Figure DI-1. The satellite's position 
is measured in terms of its original location, which is moving in the initial circular orbit. 
The coordinates (x,y) are always referenced to this initial, though moving, point. The x 
coordinate is collinear with the position vector of the initial point. The y coordinate is in 
the direction of motion of the initial point and aligned with the local horizontal. 

Figure IH-l shows the changing geometry of the scenario. Point 1 is the position 
in the initial orbit at which the maneuver begins and has unit direction vectors xi and yj. 
At some later time, the origin has moved to point 2 and the spacecraft position is 
measured in terms of X 2 and y 2 . The spacecraft has moved from the initial circular orbit 
to the final circular orbit. The angle <J> is the angle between the vector to the spacecraft 
position and the vector to the current position of the x,y origin. 
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Figure HI-1 Geometry for Hill’s Equations of Motion 


The equations of motion are developed by analysis of relative positions and 


velocities as described by Vallado. [Ref. 5] After reducing to first order differential 

equations, the two dimensional equations are given by 


x = V x 

(m-i) 

T 

V x = 2(iN y + 3co 2 x + f x = 2oN y +3co 2 x+ —sin(f) 
y J m 

on-2) 

>* 

II 

(m-3) 

T 

V y = -2 oN x + f y = -2coV x +—cos(^) 

3 y m 

(m-4) 
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-T 


(m-5) 


m = 


2. Normalization 

A normalization was performed using the initial orbit radius and velocity as 
references. A reference radius, r re f, was defined as the initial orbit radius. A reference 
time, t re f, and a reference velocity, v re f, defined by 


t = — 

l ref -l 


(m-6) 


v ref = , — 

Vnf 


(m-7) 


permit the equations to be normalized as 


x = V* 

__T 

V x = 2 + 3g> 2 x +— sin(s) 

m 

y = v y 

V =-2m x +=cos(s) 
m 


m = - 




v v «y 


(m-8) 

(m-9) 

(m-io) 

(m-ii) 

(m-i2) 


B. NONLINEAR ORBITAL EQUATIONS - BRYSON HO EXAMPLE 
1. Orbital Equations of Motion 

Bryson and Ho [Ref. 6] present a nonlinear orbit transfer formulation using 
constant thrust for a fixed time. This problem has been solved and a "known" answer is 
used as a baseline for comparison. The orbit geometry is presented in Figure III-2. The 
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thrust angle, <j>, has been replaced with £ to correlate with Hill's equations previously 
developed. 



Figure IH-2 Geometry for Nonlinear Equations of Motion (from Bryson and Ho) 


A modified formulation allowing for variable thrust results in the following 
normalized equations of motion 


r = u 

(m-i3) 

v 2 u T sin(f) 
u =- -hr + -— 

(m-i4) 

r r m 

uv T cos (s) 

V =- 4- - — 

(in-15) 

r m 

. -T 
m =- 

(m-i6) 


where 0 < T < 0.1405 (variable thrust formulation) 
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Note that in Bryson and Ho, the thrust is constant at T = 0.1405 

C. COMPARISON METHODOLOGIES 

The linear and nonlinear formulations discussed above were compared for two 

different problem formulations; free final time and fixed final time. In the free final time 
formulation, the objective is to minimize the transfer time from one circular orbit to 
another circular orbit using thrust and thrust angle as controls. In the fixed final time 
formulation, the objective is to maximize the orbit radius using constant maximum thrust 
and only the thrust angle as a control. 

The case to be studied consists of an orbit transfer from one circular orbit to 
another. The nondimensional orbit has an initial orbit radius of one (r=l) and an initial 
transverse velocity of one (v=l) in non-dimensional units. The thrust is modeled as a 
nondimensional quantity equal to 0.1405 for the constant thrust scenario. In the cases 
where thrust is a control variable, it is limited between zero and 0.1405. 

The final time conditions ensure that the final velocities V x (tf) and Vy(tf) 
correspond to a circular orbit and must be written in terms of the velocity with respect to 
the central body, such as the earth. These final time conditions are 

V x (t f ) = • sin(£) + 6> 0 y (t f ) (HI-17) 

V y 0-f ) = -J~ ' COS (£) — — ®0 X (tf ) (m-18) 

The cost function (costfii) is the final radius. Maximizing the final radius is the 
same as minimizing the negative final radius. 
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D. FIXED FINAL TIME - CONSTANT THRUST COMPARISONS 


The fixed final time comparisons are based upon the final radius obtained using 
Hill’s equations and the Bryson and Ho problem. The analysis began by running the 
same problem as Bryson and Ho, that is, continuous thrust with final time (tf) equal to 
3.32 to see if the same results were obtained. The standard answer given by Bryson and 
Ho is a final radius of 1.525. By comparison, DIDO generated a final radius of 1.52 for 
the Bryson Ho formulation and 1.4964 for Hill’s formulation. This validated the Bryson 
Ho program files. 

Further solutions were obtained for times less than the original final time of 3.32. 
The results are provided below. The percent difference is the amount by which the Hill’s 
solution differed from the Bryson Ho formulation. These answers are for n=60 LGL 
points. 


Table DI-1 Fixed Final Time Comparison from NPSOL 


Final Time 

Final Radius 

Percent 

Difference 

Bryson Ho 

Hill 

3.32 

1.52 

1.4964 

1.55 % 

2.5 

1.2772 

1.3169 

3.11 % 

2 

1.1568 

1.1755 

1.62 % 


The results are very good. It appears the linear equations provide answers that 
closely approximate the nonlinear equations. Typically, linear equations are simpler and 
faster to solve. The use of Hill’s linear equations may provide a suitable substitute for the 
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nonlinear equations to obtain quick approximations for the Bryson and Ho formulation. 
Similar performance against other problem formulations is not assumed. 

Figure HI-3 contains the results of the Bryson Ho formulation for tf=3.32 and 
Figure HI-4 contains the results of the Hill’s formulation for the same final time, tf=3.32. 

The radius and thrust angle curves are similar. The Bryson Ho formulation 
provides the radial and transverse velocity states. The Hill formulation provides the x 
and y states which are plotted together. The angle <j> is derived from the x and y states. 
The radius plot in Figure m-4 is derived from 



t t 


Figure m-3 States and Controls, tf =3.32 (Bryson Ho Formulation) 





Hills, n = 30 T1=3.32 Cost = -1.4958 



t y 


Figure HI-4 Hill’s Equations States and Controls, tf =3.32 

E. FREE FINAL TIME COMPARISONS 

In the free final time analysis, the final radius is fixed and the program optimizes 
(minimizes) the amount of time to reach that orbit. Both thrust and thrust angle are 
control variables. Unlike the constant thrust used in the fixed time analysis, the thrust 
may vaiy within the constraints described by 
0<T< 0.1405 

The results of several runs at different final radii are given in Table ni-2. 
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Table Ef-2 Free Final Time Comparison 


Final Radius 

Final Time 

Percent Difference 

Bryson Ho 

Hill 

1.05 

1.1455 

1.1922 

4.08% 

1.1 

1.6182 


2.32% 

1.3 

2.5917 

2.3829 

8.06% 

1.525 

3.3194 

3.4944 

5.28% 

2 

4.6504 



2.5 

5.8054 



3 

7.8518 



3.5 

7.5065 




A solution to the Hill's equations for final radius larger than 1.5 could not be 
obtained. The difference between the Bryson Ho solutions and the Hill's solutions 
increased slightly as the final radius increased. This indicates that the linear Hill's 
equations appear to have a limitation beyond which they are not a reliable replacement for 
the more robust nonlinear equations. A larger problem was the inability to obtain a 
solution using the Hill formulation for radii larger than 1.5. 

Examples of the states and controls for both formulation types are shown in 
Figures HI-5 and HI-6. The free final time results using the Bryson Ho formulation in 
Figure ID-5 compare very well against the fixed final time results given in Figure HI-3. 
The radius, thrust, thrust angle, and mass profiles are nearly identical. The results from 
the Hill formulation shown in Figure HI-6 show trajectories with the same general shape 
as the Bryson Ho formulation. From this analysis, it appears that Hill’s linear equations 
provide reasonable approximations to the Bryson Ho formulation for radii less than 1.5. 
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rv. DIRECT METHOD ANALYSIS 


A. ANALYSIS METHODOLOGY 

The problem formulated in Chapter II was solved numerically and then checked 

against Pontryagin's Minimum Principle (PMP). The numerical solution was obtained 
using MATLAB and NPSOL in the MATLAB environment. [Ref 11] The outputs of this 
solution are then used to determine if the optimal control solution, u*, satisfies the PMP. 

1. Numerical Computation 

The NLP resulting from the LGL pseudospectral discretization was solved using 
the NPSOL software. The NLP problem must be stated in the form: minimize f(x) 
subject to the constraints / < k(x) < u where 

k(x) = Ax 
K c(x) 

The vector x is a set of states and controls (called xopt in this analysis), f(x) is a 
nonlinear function, A is a matrix that accounts for linear constraints, and c(x) is a vector 
of nonlinear functions/constraints. The functions f(x) and c(x) are assumed to be smooth, 
i.e., at least twice-continuously differentiable. 

The problem, as defined above, was input into four different MATLAB script 
files: optmainfixed3.m, optconfixed.m, optobj.m, and optinitialfixed.m. The file 

optmainfixed3.m is run from the MATLAB command line and calls the others, along with 
the NPSOL software, when needed. 
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The optmainfixed3.m file is formatted as specified in the NPSOL user's manual. 
The entire file is given in Figure IV-1. The key parts are marked with large bold Roman 
numerals. Part I defines the A matrix. This matrix accounts for the linear constraints 
and is generally of the form A * xopt < b, but for this problem it is A* xopt = 0. The 
matrix Aisan(mxk*n) dimensional matrix where m is the number of linear constraint 
equations, n is the number of LGL points, and k is the total number of states plus controls 
(7 in this problem). For this problem, the A matrix is used for the linear periodic 
constraints: r(0)-r(tf)=0, v(0)-v(tf)=0, and y(0) - y(tf) = 0 . 

The xopt vector contains the values of each state and control at each LGL point. 
This vector will eventually contain the optimized state and control histories at the n LGL 
points (see below). 

r(l) 

r(n) 

v(l) 
v(n) 

gamma (1) 

gamma (n) 
m(l) 

m(n) 

T(l) 

T(n) 
eps(l) 

eps(n) 
theta (1) 

theta (n) 
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clear all 

global n Dn xx w t Tf 
global Ve B; 

% Define constants 
r0=300; 

vO=sqrt(398600.5/(6378.15+rO))*1000; % initial velocity 

Isp=300; 

Ve=Isp*9.8l/v0; 

Tf=112.6 

B=40900; % Generic Spacecraft, B= (2*m0)/ (r0*rho0*area*Cd) 

n=25 % Number of LGL points 

[Dn,xx,w] =dif fmat (n) ; 
t=(Tf/2)*(xx+1); 

% There are linear periodic constraints; therefore need A 

A=[1,zeros(l,n-2),-1,zeros(1,6*n); % r(l)-r(n)=0 

zeros(1,l*n),1,zeros(l,n-2),-1,zeros(1,5*n); % v(l)-v(n)=0 

zeros(1,2*n),1,zeros(l,n-2),-1,zeros(l,4*n)]; % gamma(1)-gamma(n)=0 

% Lower and upper limits for [r;v;gamma;m;T;eps;theta;A;c] 

% r(1)=1, theta(1)=0 

% r ; v ; gamma ; mass ; 

1=[1;zeros(n-1,1);zeros(n,1); -(pi)*ones(n,1);1; zeros(n-1,1);... 

zeros(n,1);-pi*ones(n,1);zeros(n,1);0;0;0;zeros(5*n, 1) ] ; 

% thrust ; eps ; theta ; A ;5 c eqns ] 

% r ; v ; gamma ; mass ; 

u=[1;inf*ones (n-1,1);inf*ones(n,1); (pi)*ones(n,1);ones(n, 1);- 

5*ones(n,1);pi*ones(n,1);0;inf*ones(n-1,1);0;0;0;zeros(5*n,1)]; 

% thrust ; eps ; theta ; A ;5 c eqns 3 

[xoptO] = optinitialfixed; 
funobj = ' optobj'; 
funcon = 1 optconfixed*; 
verifyLevel=0; 
derivativeLevel=0; 

[xopt, f,g, c, CJac, inform, lambda, iter, istate] =npsol (A, l,u,xoptO, funobj , fun 
con, verifyLevel, derivativeLevel) ; 

r=xopt(lsn); v=xopt(n+1:n*2); 

gamma=xopt(2*n+l:n*3) ; m=xopt(3*n+l:n*4); 

T=xopt(4*n+l:n*5); eps=xopt(5*n+l:n*6 ); 

theta=xopt (6*n-hl :n*7) ; _ 

Figure IV-1 optmainfixed3.m MATLAB script file 





The next significant part of optmainfixed3.m 9 part II, sets the lower and upper 
limits (1 and u) for each state, control, and constraint equation. For a free final time 
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problem, the lower and upper limits of the final time are also needed. The lower and 
upper limits on the states and controls are set to meet the constraints given in Chapter II 
and repeated here for clarity. Note that the limits at each LGL point must be specified. 
The radius initial value was set equal to one by setting both the lower and upper limit 
equal to one for r(l). The rest of the limits on r and the other states were set as wide as 
possible. For instance, the velocity has lower limits of zero and upper limits of infinity. 
Basically, no limits at all. The constraints (A and c equations) are set equal to zero by 
specifying the lower AND upper limit as zero. 

State constraints (r,v,y,m,0) 

r(0) - r(tf) = 0 (or r(0) = r(tf) = 1 for the pinned boundary condition) 
v (0) - v (tf) = 0 
y(0)-y(tf) = o 



0 < m < 1 (normalized mass) 

o<e 

Control Constraints (T,s) 

0 < T < 5 (normalized thrust) 

-ic <8 <n 

Part IQ contains the command line that actually calls the NPSOL software. Inputs 
include the A matrix, the lower and upper limits (l,u), the initial guess xoptO, the 
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objective function or cost (funobj), the constraint function (funcon), the verifylevel and 
derivativelevel. 

Several outputs are also generated. These include the optimal state and control 
histories (xopt), the final value of the objective/cost (f), an array of the objective gradient 
(g), the final values of the nonlinear constraint functions (c) and the final values of the 
Jacobian matrix of the nonlinear constraints (CJac). The value "inform" reports the result 
of the call to NPSOL. An inform^) is desirable and means that "the iterations have 
converged to a point that satisfies the optimality conditions". The term "iter" is the 
number of major iterations performed. The term "istate" describes the status of the 
constraints. [Ref. 10] 

The output "lambda" is used to generate the costates. These costates are 
important in verifying the necessary conditions of optimality. 

The five normalized state Equations (II-8 through 11-12) are the five constraint 
equations ("c" equations) contained in optconfixed.m as shown in Figure IV-2. The first 
part of optconfixed.m is a simple exponential density model. Since optconfixed.m is 
called repeatedly, the density values change as the optimal radius history vector, r, is 
found by NPSOL. The constraint ( or "c") equations are the nonlinear state equations in 
the form x - f = 0 and formatted using the differential matrix Dn. 

Each scalar constraint must be met at every point in time, therefore, it is written as 
a vector that contains the values of the constraint at each LGL point. As an example of 
the formulation, the r state equation (Equation 13-14) is formatted as a constraint, and then 
written in the required format. Note that r, v, and gamma are all nxl vectors. 
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(state equation) 


(H-14) 


r = v • sin(y) 

r - v • sin(/) = 0 (constraint equation) 

c(l :n,l)=(2/ tf)*(Dn*r)-(v.*sin(gamma)); (formatted constraint equation) 


function [c] = optconfixed(xopt) 

% NPSOL Implementation of non-linear equations of motion and constraints 
% Note, this only works for case where initial altitude is 300 km. 

% Otherwise, you must change rinitial below. 

global n Dn xx w t; 

global B Ve Tf; 

r =xopt(l:n); 

v =xopt(n+1:n*2); 

gamma =xopt(2*n+l:n*3); 

m=xopt(3*n+l:n*4); 

T=xopt(4*n+l:n*5); 
eps=xopt(5*n+l:n*6); 
theta=xopt(6*n+l:n*7); 

% Calculate density for use in drag equations 

% Changes density only when r is greater than specified 

rho=ones(n,1); % initially assume normalized density =1 for all r 

rinitial=6678.15; 

jl=find(abs(r)>(rinitial-25)/rinitial); % returns indices where true 
rho(j1)=1.87e-ll*exp(-1*(r(jl)*rinitial-6678.15)/50.3)/(1.87e-ll); 
jl=find(abs(r)>(rinitial+25)/rinitial); % returns indices where true 
rho(jl)=6.66e-12*exp(-1*(r(j1)*rinitial-6728.15)/54.8)/(1.87e-ll); 
jl=find(abs(r)>(rinitial+75)/rinitial); % returns indices where true 
rho(jl)=(2.62e-12)*exp(-1*(r(j1)*rinitial-6778.15)/58.2)/(1.87e-ll); 
D=rho.*(v. a 2); 
t=(Tf/2)*(xx+1); 

c(1:n,1)=(2/Tf)*(Dn*r)-(v.*sin(gamma)); 

c(l*n+l:n*2,1)=(2/Tf)*(Dn*v)-(-(1./r. A 2).*sin(gamma) 

+(T.*cos(eps)-D)./ (m*B)); 
c(2*n+l:n*3,1)=(2/Tf)*(Dn*gamma)-(((v. A 2)./r 

-(l./r. A 2)).*(cos(gamma))./v + (T.*sin(eps))./(m.*v*B)); 
c(3*n+l:n*4,1)=(2/Tf)*(Dn*m)-(-T./(Ve*B)); 

c(4*n+l:n*5,1) = (2/Tf)*(Dn*theta)-(v.*cos(gamma)./r) ;_ 


Figure IV-2 Nonlinear Constraints File, optconfixed .m 


The objective or cost function is incorporated into the file optobj.m. This file, as 
shown in Figure IV-3, simply creates the required function file for input into the NPSOL 
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command line in optmainfixed3.m. The input is the xopt vector and the output is the 
value of the cost (objf). 


function [objf] = optobj(xopt) 

% Created by Capt Larry Halbach, Naval Postgraduate School 
% 3 Sep 99 

global n Dn w 

r =xopt(l:n); 
v =xopt (n+1: n*2) ; 
gamma =xopt(2*n+l:n*3); 
m=xopt(3*n+l:n*4) ; 

T=xopt(4 *n+l:n*5) ; 
eps=xopt(5*n+l:n*6) ; 
theta=xopt(6*n+l:n*7) ; 

objf = sum(w.*T)/2;_ 


Figure IV-3 Cost Function File, optobj.m 


One more file was used to create the initial guess. The initial guess can impact 

the results and therefore should be as good as possible. A function file called 

optinitialfixed.m , shown in Figure IV-4, was created to generate the initial guess based on 

the results of a previous run, typically at smaller n. If a previous run was not available, 

the following initial guesses were used. 

rO=ones (n, 1) vO=ones(n, 1) 

gammaO=zeros(n, 1) m0=(1-.0001 *linspace(l, 10,n))’ 

T0=5*ones(n,l) epsO=zeros(n,l) 

theta0=(linspace(0,Tf0,n))’ Tf0=l 12.6 
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function [xoptO] = optinitialfixed 

% This function generates and returns the initial guess for optmain 
% to run NPSOL. It splines the data from a run with lower n into 
% a guess for a run with larger n. The filename in the load command 
% must be changed to the appropriate name. 

% 

% Capt Lawrence Halbach, Naval Postgraduate School 
% 13 Sep 99 

global xx n 

load j enNPS0LN20 t xopt; % ensure correct filename and old value 

old=20; % value of n from loaded filename 

to=t; % simply renames the t vector 

t= (Tf/2) * (xx+1) ; % recreates the t vector, same as in main program 

% Generate initial guesses 

rO=spline(to,xopt(l:old),t) ; 
vO=spline(to,xopt(old+1:2*old),t); 
gammaO=spline(to,xopt(2*old+l:3*old),t); 
mO=spline(to,xopt(3*old+l:4*old) , t) ; 

TO=spline(to,xopt(4*old+l:5*old),t); 
epsO=spline(to,xopt(5*old+l:6*old),t); 
thetaO=spline(to,xopt(6*old+l:7*old),t); 

clear xopt; 

% Incorporate into initial guess column vector 

xoptO= [rO; vO; gamma0;m0; TO ;epsO ; thetaO] ;__ 


Figure IV-4 Initial Guess Function, optinitial.m 


B. RESULTS 

1. Comparison to Jensen’s Results 

The state and control histories obtained using NPSOL as the NLP solver are 
presented in Figure IV-5 for n=30 and tf=112.6. The results are compared to the results 
obtained by Jensen using constr.m as the solver for the same fixed final time of tf =112.6 
[Ref. 3]. 
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Generic Satellite in NPSOL n = 30 Tf = 112.6 Cost = 0.78453 




Figure IV-5 States and Controls using NPSOL Solver (tf =112.6, n=30) 


The radius increases to a maximum of 1.0034. This equates to an increase in 
altitude of 22.7 km over the original altitude of 300 km. The radius slowly returns back 
to the beginning altitude and completes this cycle in the specified tf =112.6 time units. 
Jensen’s results also peaked at 1.0034 with a similar profile. Note that the initial radius 
(at t=0) is fixed at 1.0 as discussed previously and shown by the code in Figure IV-1. 

The velocity profiles were also similar. Both resemble an inverted radius profile 
and return to the original value at the required final time. 
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The flight path angle (y) is nearly zero throughout the entire time period. This 
indicates that the orbit remains nearly circular. Note that the units in the plot are degrees. 
In radians, the change in y is a barely discemable 0.00024. 

The mass profile is a direct result of the thrust profile. When thrusting, the mass 
rate of decrease is directly proportional to the thrust as given by the state Equation II-9. 
When not thrusting, the mass remains constant. The final mass of the spacecraft, 0.9942, 
is a 0.58 percent decrease from the mass at the beginning of the maneuver. 

The thrust profile behaves in ''bang-bang" mode. This result was not found by 
Jensen using MATLAB's constr.m file. Jensen concludes that the thrust profile is “a 
smooth, continuous throttle bum.” [Ref. 3] Under bang-bang control, the thrust will 
either be zero or maximum. In Figure IV-5, the thrust increases quickly to maximum 
thrust (T max =5), decreases very quickly to zero and remains zero until near the end. The 
final thrusting at the very end of the time period is also at Tmax- This thrusting is 
apparently used to circularize the orbit and meet the final time conditions. The impact of 
this final thrusting is particularly visible in the plot of the flight path angle. Further 
analysis on the thrust profile is conducted in Chapter VI. 

The thrust angle shows little variation during thrusting. In Figure IV-5, the thrust 
angle quickly goes to near zero as the thrust increases to T^ and then remains near zero. 
The thrust angle appears to change randomly during periods when the thrust is zero. This 
is a numerical issue since the steering is physically irrelevant when the thrust is zero. 

The cost, as given by Equation 11-20, is a ratio relating the amount of propellant 
used by the minimizing trajectory to the amount of propellant used by an FKT performed 
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at the original altitude. A cost of less than one indicates that this profile is more efficient, 
i.e. uses less propellant, than a low FKT maneuver for the same period. 

The profile in Figure IV-5 yields significant savings over a low-FKT and agrees 
well with the findings by Jensen. For this scenario, the cost equaled 0.78453. This is a 
savings of 21.5 percent over an FKT performed at the original altitude. Jensen obtained a 
cost of 0.7837 which is a difference of l/10th of one percent. 

2. Effect of Final Time Constraint 

The chosen final time greatly affects the cost index. The longer the period, the 
less the cost. Figure IV-6 shows how the costs vary with final time. This clearly shows 
the expected result that the cost, a ratio of the cost of the periodic reboost to the cost of a 
low-FKT, decreases with a longer period. 

For example, the cost of 0.39729 for the tf =700 Case is much less than for the 
shorter period discussed above. The longer period yields a higher maximum altitude 
from which it takes longer to descend. At the higher altitude, velocities are less and the 
densities are much less. The minim um density experienced at the maximum radius in 
Figure IV-5 is 0.2622 which equates to 4.9031e-12 kg/m 3 . The minimum density for the 
tf =700 case is 0.1915 which equates to 3.581 le-12 kg/m 3 for this analysis. These factors 
decrease the cost by nearly 50 percent compared to the tf =112.6 case and by over 60 
percent compared to a low altitude FKT. 
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These trajectories were then compared to a Hohmann transfer, which traditionally 
is thought to be the most efficient transfer method. As described in Chapter E, a mid- 
FKT is similar in performance to a Hohmann transfer [Ref 1]. Therefore, the cost ratio 
of the periodic reboost method to the mid-FKT is almost the same as the cost ratio to a 
Hohmann transfer. 

By dividing the periodic cost ratio plotted above by the mid-FKT cost ratio, a new 
measure of cost was obtained which is the ratio of the periodic reboost to the mid-FKT. 
Similar to the above analysis, a ratio less than one means that the periodic reboost has 
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lower costs. The following paragraph describes by example how the ratios were obtained 
for the case of n=30 and tf =112.6. 

From Figure IV-5, the middle altitude was determined to be 1.0017 which is 
halfway between the initial radius (r=l) and the highest radius (r=1.0034). The 
normalized atmospheric density at r=1.0017 is 0.7980 and the FKT velocity is 0.9983 ( = 
sqrt(l/r)). This creates a drag, D, at this altitude of 0.7953. By FKT definition, the thrust 
equals the drag. When the thrust is set equal to this value at all of the LGL points, the 
cost was determined to be 0.7953 with respect to a low-FKT as shown by 


1 f 

FKT Cost =— fTdt = 

n 



The ratio of the periodic reboost to the low-FKT was 0.78453. Therefore, the 
ratio of the periodic reboost to the mid-FKT is 0.78453/0.7953 which equals 0.9865 for tf 
= 112.6. Since this is less than one, it appears the periodic reboost is more efficient than 
the mid-FKT. Table IV-1 and Figure IV-6 contain the results from runs conducted with 
varying final time (tf). The mid-FKT cost was recomputed for each run as described 
above for tf = 112.6. As stated above, the cost of a mid-FKT is comparable to a 
Hohmann transfer. Therefore, the periodic reboost appears to be more efficient than the 
mid-FKT and the Hohmann. 
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Table IV-1 Mid-FKT Cost Analysis 


Final Time 

Periodic Reboost 

Mid-FKT 

Periodic Reboost 

Low-FKT 

Low-FKT 

Mid-FKT 

112.6 

0.78453 

0.7948 

0.987 

130 

0.7688 

0.7843 

0.980 

150 

0.72328 

0.7396 

0.978 

200 

0.66194 

0.6861 

0.965 

300 

0.57618 

0.6148 

0.937 

400 

0.51336 

0.5177 

0.992 

500 

0.46492 

0.4757 

0.977 

600 

0.42718 

0.4412 

0.968 

700 

0.3954 

0.4174 

0.947 


Analysis for periods less than 112.6 yielded suspect profiles in which the state and 
control histories are not entirely reasonable. As an example, the states and controls for tf 
=80 show erratic behavior as shown in Figure IV-7. The radius and velocity histories 
have unrealistic spikes near the beginning and the end. Therefore, these cases are not 
included in the analysis. 

Interestingly, the plot of the mid-FKT cost comparison in Figure IV-6 shows a 
local minimum near a final time of 300. A time free analysis was performed to examine 
this possibility. 

The free final time run had the final time limited to between 100 and 400. A local 
minimum would result in a final time less than 400. In this case, the final time increased 
to the maximum limit of 400. The cost of 0.51398 is nearly identical to the cost 
(0.51336) obtained from the fixed final time run of tf =400. The mid-FKT/low-FKT cost 
ratio, calculated as above, equaled 0.5191 which yields a cost index of 0.990. The free 
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final time results agree well with the fixed final time results for tf=400. Based on this 
analysis, a local minimum at 300 is unlikely. 


Generic Satellite in NPSOL n = 30 



Tf= 80 Cost = 0.8312 




Figure TV-7 State and Control History, tf =80 


3. Effect of Radius Constraint 

Many spacecraft are operated within an operational orbital belt that places both an 
upper and lower altitude restriction on the orbit This prevents complete implementation 
of the optimal periodic reboost. Consider the case of tf =700. The unconstrained 
maximum altitude is 1.0123. Setting the upper radius limit to 1.008 (equivalent to 53.4 
km) ensured a bounded constraint. The results are presented in Figure IV-8. 

The profile begins the same way as the unconstrained case with maximum thrust 
and rapid increase in radius. But, when the upper limit is reached, the thrust decreases 
and the radius flattens out until much later when it begins the familiar decay due to drag 
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back to the original altitude to meet the final time requirement. While maintaining the 
maximum allowed altitude, the spacecraft is flying an FKT trajectory, since at the 


maximum altitude, r=1.008, and v = 



0.99602. Thrust stabilizes at normalized 


thrust equal to normalized drag as given by 

D = pv 2 = 0.33457 * .99602 2 =.3319 

At the original altitude, an FKT trajectory would require T=l. But at the higher 
altitude with its lower velocity and lower density, the required thrust is much less. This 
profile results in a higher cost than the unlimited altitude case: 0.44157 vice 0.39729. 


Generic Satellite in NPSOL n = 50 
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Figure IV-8 States and Controls for Altitude Limited Case (tf =700, n=50) 
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4. 


Costates and Necessary Conditions 


In this section, the outputs from the direct method analysis are checked against the 
necessary optimality conditions. 

The necessary conditions for optimality, as described in Chapter II are 


„ <T dH 

Costates: A =- 

dx 

(D-l) 

„ . SH A 

Optimality: — = 0 

du 

(D-2) 

o 

A1 

w « 

1 

(D-3) 

T ,, x 0M ( 

Transversahty I: A(tf)= + V f 

oXf {dx f J 

(n-4) 

Transversality II: (if final time free) 


H( k ) + 8M(tf) + vJ d l (tf Ko 

otf ctf 

(D-5) 

where, 


H = the augmented Hamiltonian = L + A T f + p. T g 



in which x = f from the state equations and g is the control inequality vector 


given by 


Si' 


T - Tj^ 


0 ' 
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-T 
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► < < 
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S3 


e-n 
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.84, 


-e-n 

V J 
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X 

fi g has the relation that p > 0 when g = 0 
and ji= 0 wheng<0 

Since g 3 and g 4 are always less than zero, then p. 3 = 0 and ji 4 = 0. 


Therefore, ji t = i ^ 
[Ml. 


The cost was given in Chapter II as 


J = — (* f Tdt 
t f •*> 


Also, X = column vector of costates = < 




^tf) s final conditions. 


u = vector of controls = ■ 


and x = vector of states = \ 


The first part of the MATLAB script file optfixresults.m, takes the NPSOL output 
"lambda" and generates the costates related to the state equations by using the equation 
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(IV-1) 


X(ti) 


%LP 

Wj 


where Wi are the weights at the LGL points. [Ref. 4] 

The vector lambda from NPSOL contains additional values, besides the costates. 
The first 7*n values (n=number of LGL points) contain the lambdas associated with the n 
initial conditions on each state (r,v,y,m,0) and each control (T,s) in the following order: 
(r,v,y,m,T,s,0). The next three components, lambda(7*n+l) to lambda(7*n+3), 
correspond to the linear constraints. The next 5*n elements contain the costates for the 
five state equations for (r,v,y,m,0). These are the values of interest for analysis of the 
necessary conditions given above. 

The code to accomplish this in the MATLAB file is given in Figure IV-9. 


lamr=lambda(7*n+4:8*n+3)./(w*(Tf/2)); 
lamv=lambda(8*n+4:9*n+3)./(w*(Tf/2)); 
lamgamma=lambda(9*n+4:10*n+3) ./(w* (Tf/2)) ; 
lammass=lattvbda (10*n+4 : ll*n+3) . / (w* (Tf/2)) ; 

lamtheta=lambda(ll*n+4:12*n+3),/(w*(Tf/2)) ; _ 

Figure IV-9 Costates Generation Code 

These costates were then plotted against time. Figure IV-10 shows the plots for 
the results of a run at tf=l 12.6 and n=30. The costates obtained from this analysis were 
noisy. A digital filter designed by T. P. Thorvaldsen and based on the research of Ross, 
Fahroo, and Thorvaldsen was applied to smoothen the costates. The filtered costates are 
presented in Figure IV-11 and are used in the following analysis. The filtered costates 
contrast sharply with the unfiltered costates. The large cyclic changes in magnitude are 
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Figure IV-10 Plots of Costates 


46 





Costates, n = 30 Tf= 112.6 


Cost =0.78453 
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Figure IV-11 Filtered Costates for n=30 


An analysis was performed using these filtered costate estimates. The analysis 
begins with the Hamiltonian, which is given as: 


H = Hamiltonian = L +X T f + ja T g 

gl 

.§2 


H= — + W 
t f 




(IV-2) 
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TT T „ n f-sinv Tcos^-D | 0 v 1 cos^ Tsin e 

H = —+Vsmr + /lv —H-+ +4-^ L + —— 


+ X 0 ?^ +/I (j- T ) +fl2 (-T) 

r 

Note: All of the values used in these equations are the normalized values. The 
bar over the symbol is not included for clarity of presentation. 

The PMP states that the Hamiltonian, H, is a minimum at the optimal control u*. 
This applies at every point in the trajectory. 

The optimality conditions for the controls are given by 


*=0 

du 


(IV-3) 


®. r * ®i.[j o] 

du ST 5s 


(IV-4) 


SH 1 , cose sins , f 1 A 

— = — + Ay —— + 'V —~ _ —~ + Ml ~ 7*2 - 0 

ST tf mB mvB 1 v P B 1 


(IV-5) 


5H . Tsrnf . Tcosf 

— = -Z v -+ 7L-= 0 

ds mB mvB 


(IV-6) 


In —, pi equals zero except where the thrust, T, equals Tmax and p 2 equals zero 
ST 

everywhere except where T=0. The thrust should change from Tmax to T=0 or vice versa 
whenever the switching function, S, changes sign, where S is given by 


0 1 , coss sins [ 1 

tf mB r mvB v„B 


(IV-7) 


To check this condition, the switching function and thrust profile, scaled to 


171000 th , are plotted together against time. In Figure IV-12, the original thrust profile 
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from Figure IV-5 is plotted with the switching function that results from the unfiltered 
costates. As can be seen, as the switching function crosses from negative to positive near 
time 15, the thrust drops to zero. Later, the switching function becomes negative for a 
moment and the thrusting increases (to Tmax). It appears the solution obeys the switching 
function given by Equation TV-7. 


x 10‘ 3 



Figure IV-12 Switching Condition, S, and Thrust/1000 

In Figure IV-13, the filtered costates are used to generate the switching function. 
Based on this switching function, the thrust profile was determined and is also plotted in 
Figure IV-13. The filtered costates yield a switching function that leads to a bang-bang 
thrust profile. The thrust is a maximum when the switch is negative and goes to zero 
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when the switch becomes positive. The switching function is discussed further in 
Chapter VI. 

-3 

x 10 



Figure IV-13 Switch and Thrust Using Filtered Costates 

dH 

Figure IV-14 contains the graph of the other optimality condition, -. 

os 

According to Equation IV-6, this should equal zero at all times. The results generally 
agree with this requirement except for large errors at the beginning and final times. The 
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cause of these errors is unknown, but the deviations occur in every plot of the costates 
and necessary conditions. Overall, the optimality conditions appear to be met. 


x io' 3 



Figure IV-14 Optimality Condition, — = 0 

ds 


The second partial of the Hamiltonian with respect to the unconstrained controls. 


0 2 H 

dn 2 


> 0, is also an important tenet of PMP. Since thrust is constrained, only the second 


partial with respect to the thrust angle is analyzed and is given by 
3 2 H _ , Tcosf , Tsins 

-— — —/i v - A> 

ds mB mvB 


(IV-8) 


51 





Figure IV-15 provides a plot of 


ds 2 


It appears that the steering angle does 


d 2 U 


meet the requirement that —— > 0 


ds 2 



Figure IV-15 Plot of ^ 

ds 


The Hamiltonian, H, should be zero for optimality. Figure IV-16 contains a plot 
ofH as given by 


TT T , , f-siny Tcoss-D^l , 

H = — + A r vsm^ + A v —-- \+A,y 


f 2 

V 


cos y T sin s 
- L.+ - 


Xm T v cc 


■rmrTi 






Figure IV-16 Plot of Hamiltonian 


As shown in Figure IV-16, the Hamiltonian is nearly constant at H=0 except at the 
beginning and the end. 

The necessary conditions for the costates are given by Equation II-1 which can be 
rewritten as 


X J 


+ 


m 

dx 


= 0 


CIV-9) 


The individual components become 







In order to evaluate these expressions, the values of X must first be computed. 
This was done using the differential matrix introduced earlier. The formulation is given 
by 

A =f-l(D„*A) (W-15) 

Uf J 

The costate necessary conditions were then plotted in Figure IV-17 based on 
Equations IV-10 to IV-14. The plots should equal zero if the results match the theory. In 
general, these plots do meet the requirements fairly well except for ^ and at the 
beginning and end of each plot. 
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time 


• t 3H 

Figure IV-17 Costate Difference Plots, A H-= 0 


The transversality condition is a function of any final time cost elements, M(t f ), 


and the final time conditions, ^(tf). In this problem, the cost function is an integral 


function and does not have a final time cost element, so M(t f ) = 0. Therefore 


■ = 0, which results in the transversality condition 
dx 


^ T / t \ dM(t f ) dv|/(t f ) _^ T ay(t f ) 

f dx dx dx 


In this problem, the final conditions are as given by 
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(IV-17) 


r(t f )-r(0) 


V'(tf) = iv(tf)-v(0)^ 


[y(t f )-y(0)J 


The partial derivative with respect to the states is given by 


dtKtf) 

dx 


a(r(tf)-l) a(r(tf)-l) 

dr dv 

a(y(tf)-i) a(v(t f )-i) 

dr dv 

a(y(tf)-KQ)) a(r(tf)-y(Q)) 

dr dv 


3(r(t f )-l) 

dy 

a(v(t f )-i) 

dy 

8(r(t f )-y(0)) 

dy 


0 0 
0 0 
0 0 


(IV-18) 


which becomes: 


dx 


1 0 0 0 0 
0 10 0 0 
0 0 10 0 


The resulting transversality equation is then 


A T (t f ) 4 r 


Uy 



1 0 0 0 0 
0 10 0 0 
0 0 10 0 


A T (tf)=[y r u v v y 0 o] 

or separated into components, the transversality condition is: 

M*f) = °r 

K(t f ) = v v 
\it f ) = V r 


(IV-19) 
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C. NUMERICAL ANALYSIS ISSUES 

1. Effects of Problem Formulation 

The formatting in the "main.m" file, such as in optmainfixed3.m, impacted the 
solutions. Several cases were run in which the linear constraint matrix, A, and the lower 
and upper bounds, 1 and u, were defined differently than described above. Instead of 
including the r(l)-r(n)=0 and v(l)-v(n)=0 constraints into the A matrix, their lower and 
upper bounds of r(l), r(n), v(l), and v(n) were set equal to 1 in both the lower and upper 
bounds. This forced these to be equal, thus forcing r(l)=r(n)=l and v(l) =v(n)=l. This 
formulation differs from the original formulation in which r(l)=r(n)=l but the velocity 
terms had a lower bound of zero and no upper bound. (See part II of Figure IV-1) 
Throughout this chapter, the radius plots have initial and final normalized radii that are 
forced to equal one (r(l)=r(n)=l). 

The new formulation changed the results. The case of n=30, tf =112.6 is given in 
Figure IV-18 to demonstrate. The difference in cost is immediately noticeable. While the 
originally formatted problem gave a cost of 0.78453, this formulation yielded a higher 
cost of0.79767. 

Since cost is a function of thrust, the thrust profiles were overlaid in Figure IV-19 
to look for differences. The thrust profile labeled as "original" is the thrust profile 
obtained from the proper formulation already presented during the results discussion. 
The profile labeled as "modified" represents the formulation being discussed here. The 
profiles are generally alike, but differences do exist. 
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The initial large thrusting profiles are similar but shifted. The final thrusting 
sequence is also different. The final thrust in the modified profile spikes to T,nax=2.1 vice 
5 for the original profile. However, they both appear to be reacting according to the 
switching function. 

The settings of the lower and upper limits also have an effect. In one case, the 
lower and upper limits of the radius and velocity were changed to reduce the range of 
values to more physically realistic limits as shown below. 

Lower limits: r: ones(n,l) (r(t=0)=l, all others > 1) 

v: l;.5*ones(n-l,l) (v(t=0)=l, all others > 0.5) 

Upper limits: r: l;1.5*ones(n-l,l) (r(t=0)=l, all others < 1.5) 

v: ones(n,l) (v(t=0)=l, all others < 1) 

While these limits appeared to make physical sense and agreed with the previous 
outputs, they effected the results and resulted in costs and solutions similar to the 
"modified" formulation described above. 

2. Stability of Solutions 

The stability of the solutions was checked in two ways. First, the outputs of the 
n=30, tf =112.6 were used as the initial guess for another run. If a solution is indeed the 
optimal solution, then, when used as the initial guess, the output should be the same as 
the guess. In this case, good stability was demonstrated by the return of a nearly identical 
solution. 
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The program was also run using a slightly modified initial guess. The different 
initial guess consisted of guessing r0=ones(n,l) instead of using the radius results of a 
previous run. The results also compared well against the previous run. 

3. Effect of Increasing the Number of LGL Points 

The number of LGL points effected the solutions and sometimes did not yield any 
solution. For example, the case of tf =112.6 was run with numerous different numbers of 
LGL points. The results presented earlier resulted from n=30. Runs of n=40 and n=50 
were also performed. The radius and thrust profiles of each (n=30, 40, 50) are presented 
in Figure IV-20. The n=40 radius profile is similar to n=30 but the thrust profile shows a 
notch near t=10 and an additional thrust spike near t=105. The cost was 0.78106 which is 
0.44 percent less than for n=30. The n=50 run yielded poor state and control profiles. 
Numerous runs were done using different initial guesses but the results always showed 
the choppiness seen in Figure IV-20. The cost was 0.77395 which is 1.3 percent less than 
for n=30. 

A similar problem occurred for cases of tf =80. The radius and thrust profiles for 
cases of n=25, n=35, and n=45 are presented in Figure IV-21. As the number of LGL 
points increases the radius and thrust profiles become worse. The radius profile has 
discontinuities and the thrust profile becomes more erratic. 

The number of points also greatly affected the solution time. This analysis was 
conducted on a Sun Sparc 10 workstation using MATLAB 5.2 and NPSOL. Typical run 
times were 30 minutes for n=20, 1.5 hours for n=30, 3 hours for n=40, and 5 hours for 
n=50. 
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V. 


ANALYSIS BY AN INDIRECT METHOD 


Indirect methods of optimization rely on the necessary optimality conditions 
derived from the minimum principle. These necessary conditions are solved for the 
optimal trajectory by solving a nonlinear two point boundary value problem to obtain the 
Lagrange multipliers and costates. [Ref. 11] There are several difficulties with indirect 
methods. Firstly, the method requires that the necessary conditions be derived 
analytically. Secondly, the indirect methods require a very good initial guess [Ref. 9]. 

The output of the direct method, including the costates estimates, are used as the 
initial guess for a multiple shooting algorithm developed at NPS. This code contains 
three primary steps. First, it loads the states and costates from a data file for use as the 
initial guess. This data file was either the output obtained from the direct method or the 
output from a previous run with this software. The software then sets up the initial 
Jacobian matrix. Lastly, it uses a quasi-Newtonian iteration to get the solution to the two 
point boundary value problem (TPBVP). The user controls the iteration accuracy 
requirement and, more importantly, the number of iterations. The number of iterations 
must be large enough to reach convergence. 

Unfortunately, a converged solution was not obtained. The indirect method is 
extremely sensitive to the initial guess. Both the filtered and the original “noisy” costate 
estimates obtained from the direct method were used as the initial guesses. Neither 
resulted in a converged solution and both are probably too poor for convergence. 
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Alternatively, there might be errors in coding as no analysis was done to validate the 
process. Attempts with slightly modified costates improved convergence, but still did not 
provide a converged solution. Therefore, this chapter presents the method of analysis and 
shows the “best” solution obtained using the original "noisy" costates. These results are 
neither optimal nor feasible and, therefore, should not be compared against the results of 
the direct analysis. 

The states and costates for the initial guess are loaded from the data file, but not 
the controls. The controls must be obtained in terms of the states and costates. Figure 
V-l contains the code that calculated the thrust angle, e, and thrust magnitude, T, based 
on given states and costates. 


eps=atan (lamdg/ (v*lamdv) ) ; 

Switch=l/Tf+lamdv*cos(eps)/(ma*B)+lamdg*sin(eps)/(ma*v*B)-lamdm/(Ve*B); 

if Switch>=0 
T=0 ; 

elseif SwitchcO 
T=Tmax; 

end _ 


Figure V-l Thrust and Thrust Angle Equations 


These equations are based upon the optimality condition and switching function 
developed in Chapter IV and repeated here. 

(IV-6) 


3H , Tsins Tcosf _ A 

-— — /i-b A — U 

8s mB mvB 


1 cos s sin^ 

t f mB 7 mvB 


< 1 N 


V v e B y 


(IV-7) 
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For bang-bang control, the thrust is a maximum when the switching function is 
less than zero and a minimum when greater than zero. The switching function is 
discussed further in the Chapter VI. 

Converged solutions were not obtained. Convergence was measured by the 
Euclidean norm of the difference between consecutive solutions during the iteration 
process. A converged solution returns a value of approximately 10' 5 to 10* 6 . Values on 
the order of 10‘ 2 were obtained. The non-convergence is measurable by the norm and 
visible in the states and costates. 

The current non-converged results are presented for the fixed final time case of tf 
= 112.6. Figures V-2 and V-3 provide the best solutions obtained with a convergence 
measure of approximately 5.75*10‘ 2 . The states have an oscillatory motion. This 
oscillation decreases in amplitude as the convergence measure decreases. 

The costates of periodic states must also be periodic. Periodic states or costates 
were not obtained. The costates of states without a final time constraint should be zero 
(see Equation IV-19). In Figure V-3, the theta costate meets this requirement. The mass 
costate is close to zero but the other costates are not periodic. Note that, since a 
converged solution was not obtained, these comparisons against theory are simply to 
illustrate how non-optimal solutions may be detected visually from the costate plots. 


65 






0.0251 



•0.01 [-*---*-1 

0 30 60 80 120 


Figure V-2 Non-Convergent States an( 



Figure V-3 Non-Convergent 


mass costates 



Convergent Measure = 0.0575 















As stated earlier, the indirect solutions rely upon the optimality conditions. These 
conditions were checked against the results provided above. Figure V-4 provides the 

9H 

plot of the optimality condition on 8, — = 0. The optimality condition appears to be 

ds 

met. This is surprising since the plot is based upon data from a non-convergent result. 

The indirect method was very sensitive to the initial guess and the number of 
iterations was also critical. Obviously, the number of iterations must be enough to reach 
convergence. Low number of iterations may lead to wildly varying results. The estimate 
of the optimal solution from the direct method allowed these results to be easily 
discounted. Figure V-5 provides examples of the results obtained from very low numbers 
of iterations. Rim times were approximately one to three minutes per iteration on a 
Pentium II at 200 MHz. 
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VL ANALYSIS OF THE SWITCHING STRUCTURE 

Analysis of the switching function and the effect of varying the thrust control 
authority provided interesting results related to singular thrust arcs. [Ref. 2] Singular 
thrust arcs result when the switching function remains at zero instead of simply crossing 
it (i.e., switching). This was observed in the thrust profiles and switching functions 
obtained from the direct method. 

The effect of increasing levels of thrust was studied using DIDO. It was 
suspected that a maximum thrust limit for bang-bang control existed. [Ref. 3] Figure VI- 
1 contains the thrust plots resulting from four different thrust levels (5, 15, 50,100) for a 
fixed final time of 700 and n=50 LGL points. The cases of Tmax ~ 5 and Tmax = 15 both 
reach their maximum thrust levels. The cases of Tmax = 50 and Tmax = 100 each peak at a 
thrust level of approximately 34.5. 

Thrusting of this type is predicted based upon the totality of extremal thrust-arcs 
described by: [Ref. 2] 


T = 


0 


S<0 


T s whenever S = 0 
T S>0 


(VI-1) 


where S is the switching function and T s is the singular thrust. This highlighted the need 
to further analyze the switching function. The values of the switching function and the 
thrust profiles are plotted for each of the cases in Figures VI-2 through VI-5. Included on 
the plots are the relevant values of n, time, S, and thrust. 
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Figure VI-1 Thrust Profiles for Increasing Tmax at final time of 700 


Figures VI-2 and VI-3 contain the thrust profiles and switching function of the 
Tmax = 5 and Tmax = 15 cases, respectively. In these cases, both have periods during which 
the switching function is approximately zero at which time the thrust is neither Tmax nor 
zero, but at some varying singular thrust level, T s . The singular thrusting occurs when the 
magnitude of the switching function is less than 10‘ 6 . 
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Figure VI-2 Switching and Thrust for Tmax = 5 and tf= 700 
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The cases of Tma X = 50 and Tma* = 100, shown in Figures VI-4 and VI-5, are 
different. Each case experiences a period in which the switching function is equal to 
zero. The thrust level varies during this period and at no time does it reach Tmax* 
Though it appears to spike similar to a bang-bang type profile, it is not bang-bang. This 
is a different thrust regime known as singular thrust-arcs. [Ref. 2]. Values of interest are 
included on the plots. 
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Figure VI-4 Switching and Thrust for Tmax = 50 and tf= 700 
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Figure VI-5 Switching and Thrust for Tmax =100 and tf= 700 
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VII. CONCLUSIONS 


The direct Legendre pseudospectral method provides solutions for the periodic 
orbital maintenance problem that appear to be more efficient than a Hohmann transfer. 
Numerical analysis indicates savings of 1 to 5 percent over an impulsive Hohmann 
transfer. Savings increase when compared to a finite-bum Hohmann transfer by as much 
as 6 percent [Ref. 3]. Further analysis into the cost ratio of the periodic control to the 
mid-FKT is needed to obtain the true form of the cost plot (see Figure IV-6). 

The costate estimates resulting from the direct method were apparently noisy, 
particularly near the beginning and ending times of the period. The "noise" was filtered 
and the resultant costates compared against the optimization theory. The costates did not 
violate the optimization theory as demonstrated in Chapter IV. 

The costate estimates are apparently unreliable for this problem. This is curious 
as costate estimates for the orbit transfer problem [Ref. 4] have been obtained from 
NPSOL. As discussed in Chapter IV, the solutions are sensitive to the implementation of 
the problem. 

Analysis by the indirect method was also attempted but without success. 
Convergence of the indirect solution is difficult to obtain due to the sensitivity to the 
initial guess. The unsmooth costates obtained from the direct method solutions may be 
partially responsible for limiting convergence. 
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The thrust profiles were bang-bang with periods of thrust modulation. These 
thrust-modulated arcs appear to be singular since the switching function was equal to 
zero. Singular thrust arcs may be particularly beneficial to large flexible structures as 
thrust modulation will probably excite fewer high frequency vibrations than a bang-bang 
controller. This has potential benefits to spacecraft payloads, structures, experiments, and 
inhabitants. 

More work in this area is certainly needed. The optimal control problem is still 
not fully solved, but the solution presented here appears to provide advantages to mission 
planners. Employing the optimal orbit maintenance technique along with properly sized 
thrusters can significantly increase mission duration and reduce cost. 
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